$title  build benchmark social accounting matrix



* ------------------------------------------------------------------------------
* load implan data
* ------------------------------------------------------------------------------

* basic sets from implan
sets
  m(*)                  implan factors
  a(*)                  accounts
  inst(*)               institutions
  s(*)                  sectors
  r(*)                  regions
  f                     fixed factors                        / l, k /;

* load the IMPLAN sets
$gdxin 'build\data\merged.gdx'
$load m<f.dim2 a<t.dim2 inst<i.dim2 s<g.dim2 r=merged_set_1

* sets defining the implan factors in the sage factors
sets
  l(m)                  labor                                / empl /
  k(m)                  capital                              / prop, othp/;

* subsets for households and public institutions
sets
  h(inst)               households
  pub(inst)             public (government) institutions
  corp(inst)            corporations
  trd                   trading markets;

* stock based sectors subsumed in investment and row adjustment sector
sets
  scrap(s)              stock based sectors                  / 527, 528 /
  rowadj(s)             rest of world adjustment             / 529 /;

* load the subsets
$load h<h.dim2 pub<pub.dim2 trd<simport.dim2 corp<corp.dim2

alias (s,ss,sss), (h,hh), (inst,ii), (r,rr), (m,mm), (trd,ttrd);

* parameters from read_implan_state
parameter
  make(r,s,a,ss)        domestic industry make matrix
  use(r,s,a,ss)         domestic industry use matrix
  iuse(r,s,a,inst)      domestic institutional use matrix
  fd(r,m,a,s)           factor input matrix : industry use of factors
  fexprt(r,m,a,trd)     factor exports
  imake(r,inst,a,s)     domestic institutional make matrix
  fs(r,inst,a,m)        factor disbursement matrix
  trnsfer(r,inst,a,ii)  inter-institutional transfers
  fimprt(r,trd,a,m)     factor imports
  trnshp(r,trd,a,ttrd)  transhipments
  sexport(r,s,a,trd)    exports by sector
  iexport(r,inst,s,trd) exports by institution
  simport(r,trd,s,ss)   imports by sector
  iimport(r,trd,s,inst) imports by institution;

$load make use iuse fd fexprt imake fs trnsfer fimprt trnshp sexport iexport
$load simport iimport



* ------------------------------------------------------------------------------
* specify parameters defining the sage sam
* ------------------------------------------------------------------------------

parameters
  y0(r,s)               benchmark - output
  ld0(r,s)              benchmark - labor demand
  kd0(r,s)              benchmark - capital demand
  id0(r,s,ss)           benchmark - intermediate demand
  d0(r,s)               benchmark - domestic supply
  m0(r,s,trd)           benchmark - imports
  x0(r,s,trd)           benchmark - exports
  cd0(r,s,h)            benchmark - consumer demand
  i0(r,s)               benchmark - investment demand
  invh0(r,h)            benchmark - household investment
  g0(r,s)               benchmark - government demand
  le0(r,h)              benchmark - household labor supply
  re0(r,h)              benchmark - capital rental income
  bopdef0(r,h)          benchmark - international balance of payments
  incadj0(r,h)          benchmark - government to household lump sum transfer
  n0(r,h)               benchmark - households
  ty0(r,s)              tax rate - production
  tl0(r)                tax rate - labor use
  tk0(r)                tax rate - capital use
  tc0(r)                tax rate - consumption;



* ------------------------------------------------------------------------------
* load and balance trade data
* ------------------------------------------------------------------------------

* implan data has a small imbalance in domestic trade so the international
* trade balance is adjusted to absorb any imbalance in domestic trade

* benchmark - imports
m0(r,s,trd) = sum(ss, simport(r,trd,s,ss))+sum(inst, iimport(r,trd,s,inst));

* benchmark - exports
$onuni
x0(r,s,trd) = sum(ss, sexport(r,ss,s,trd))
              +sum(inst, iexport(r,inst,s,trd)$(not scrap(s) and not rowadj(s)));
$offuni

* parameters for trade specification and balance
parameters
  dtrd_ratio(s)         ratio of domestic imports to exports;

* ratio of domestic imports to exports (set to -1 if no domestic exports)
dtrd_ratio(s)$sum(r, x0(r,s,"dtrd")) =  sum(r, m0(r,s,"dtrd"))
                                       /sum(r, x0(r,s,"dtrd"));
dtrd_ratio(s)$(not sum(r, x0(r,s,"dtrd"))) = -1;

* if domestic exports exceed domestic imports then reduce domestic exports
* proportionally and increase international exports by the same amount
x0(r,s,"dtrd")$(dtrd_ratio(s)>0 and dtrd_ratio(s)<1) =  dtrd_ratio(s)
                                                       *x0(r,s,"dtrd");
x0(r,s,"ftrd")$(dtrd_ratio(s)>0 and dtrd_ratio(s)<1) =  x0(r,s,"ftrd")
                                                       +x0(r,s,"dtrd")
                                                        *(1/dtrd_ratio(s)-1);

* if domestic imports exceed domestic exports then reduce domestic imports
* proportionally and increase international imports by the same amount
m0(r,s,"dtrd")$(dtrd_ratio(s) gt 1) = m0(r,s,"dtrd")/dtrd_ratio(s);
m0(r,s,"ftrd")$(dtrd_ratio(s) gt 1) = m0(r,s,"ftrd")+m0(r,s,"dtrd")
                                                     *(dtrd_ratio(s)-1);

* if domestic exports are zero and domestic imports are not then zero out
* domestic imports and add the value to international imports
m0(r,s,"ftrd")$(dtrd_ratio(s) eq -1) = m0(r,s,"ftrd")+m0(r,s,"dtrd");
m0(r,s,"dtrd")$(dtrd_ratio(s) eq -1) = 0;

* if domestic imports are zero and domestic exports are not then zero out
* domestic exports and add the value to international exports
x0(r,s,"ftrd")$(dtrd_ratio(s) eq 0) = x0(r,s,"ftrd")+x0(r,s,"dtrd");
x0(r,s,"dtrd")$(dtrd_ratio(s) eq 0) = 0;



* ------------------------------------------------------------------------------
* load reamaining benchmark data
* ------------------------------------------------------------------------------

* much of the mapping from implan to sage is straight forward, adjustments are
* in the handling of institutional make and exports. government make and exports
* are assigned to sectoral production with intermediate inputs in proportion to
* government use. inventory adjustments defined as make from stock are added
* to sectoral production as capital additions, while scrap and second hand goods
* are implicitly netted out of investment. the rest of world adjustment in
* household exports is added to the production sectors and shared out among
* goods based on the final good consumption shares.

* intermediate parameters used in mapping the implan data to sage parameters
parameters
  g_make(r,s)           government production
  g_use(r,ss,s)         intermediate inputs to government production
  row_adj(r,s,h)        rest of world adjustment
  c0(r,h)               aggregate consumption
  gr(r)                 aggregate government demand;

* benchmark - intermediate demand
id0(r,ss,s) = sum(a, use(r,ss,a,s))+sum(trd, simport(r,trd,ss,s));

* benchmark - labor demand
ld0(r,s) = sum((l,a), fd(r,l,a,s));

* benchmark - capital demand
kd0(r,s) = sum((k,a), fd(r,k,a,s))
           +sum((inst,a)$(not pub(inst)), imake(r,inst,a,s)$(not scrap(s)))
           +sum((inst,trd)$(not pub(inst)), iexport(r,inst,s,trd)$(not scrap(s)
                                                                   and not rowadj(s)));

* benchmark - consumer demand
cd0(r,s,h) = sum(a, iuse(r,s,a,h))+sum(trd, iimport(r,trd,s,h));

* aggregate consumption
c0(r,h) = sum(s, cd0(r,s,h));

* rest of world adjustment - shared across households based on their share of
* regional consumption
row_adj(r,s,h)$cd0(r,s,h) = iexport(r,h,"529","ftrd")*cd0(r,s,h)/c0(r,h);

* remove the rest of world adjustment from consumption
cd0(r,s,h) = cd0(r,s,h)-row_adj(r,s,h);

* share rest of world adjustment out across commodities
x0(r,s,"ftrd") = x0(r,s,"ftrd")+sum(h, row_adj(r,s,h));

* benchmark - government demand
g0(r,s) = sum((a,pub), iuse(r,s,a,pub))+sum((trd,pub), iimport(r,trd,s,pub));

* government production
g_make(r,s) = sum((pub,a), imake(r,pub,a,s))
              +sum((pub,trd), iexport(r,pub,s,trd));

* aggregate government demand
gr(r) = sum(s, g0(r,s));

* intermediate inputs to government production
* assume that all government production uses inputs in proportion to the overall
* government consumption shares
g_use(r,ss,s)$(g_make(r,s) and g0(r,ss)) = g_make(r,s)*g0(r,ss)/gr(r);

* remove the inputs for government production from government final consumption
g0(r,s) = g0(r,s)-sum(ss, g_use(r,s,ss));

* add the inputs for government production to intermediate inputs
id0(r,ss,s) = id0(r,ss,s)+g_use(r,ss,s);

* calculate output at market prices
y0(r,s) = sum(ss, id0(r,ss,s))+ld0(r,s)+kd0(r,s)+sum(a, fd(r,"btax",a,s));

* tax rate - production -> average rate
ty0(r,s)$y0(r,s) = sum(a, fd(r,"btax",a,s))/y0(r,s);

* benchmark - investment demand
i0(r,s) = y0(r,s)+sum(trd, m0(r,s,trd))-sum(trd, x0(r,s,trd))-sum(h, cd0(r,s,h))
          -g0(r,s)-sum(ss, id0(r,s,ss));

* benchmark - household investment
* ensure balance by scaling investment by household's share of transfer to
* investment account
invh0(r,h)$sum(hh, max(0, sum(a,trnsfer(r,"inv",a,hh)-trnsfer(r,hh,a,"inv"))))
                = max(0, sum(a, trnsfer(r,"inv",a,h)-trnsfer(r,h,a,"inv")))/
                  sum(hh, max(0, sum(a,  trnsfer(r,"inv",a,hh)
                                        -trnsfer(r,hh,a,"inv"))))
                  *sum(s, i0(r,s));

* benchmark - household labor supply
* ensure balance by scaling factor demand by factor supply by household
le0(r,h) = sum((a,l), fs(r,h,a,l))/sum((hh,a,l), fs(r,hh,a,l))*sum(s, ld0(r,s));

* benchmark - capital rental income
* ensure balance by scaling factor demand by factor supply by household
re0(r,h) = sum((a,k), fs(r,h,a,k))/sum((hh,a,k), fs(r,hh,a,k))*sum(s, kd0(r,s));

* benchmark - international balance of payments
* balance of payments per household scaled by consumption share
bopdef0(r,h)  = sum(s, m0(r,s,"ftrd")-x0(r,s,"ftrd"))*sum(s, cd0(r,s,h))
                                                      /sum((s,hh), cd0(r,s,hh));



* ------------------------------------------------------------------------------
* set the year of the benchmark
* ------------------------------------------------------------------------------

sets
  t                     time           / 2016 /;



* ------------------------------------------------------------------------------
* add taxes to benchmark dataset
* ------------------------------------------------------------------------------

* labor tax rate - based on FICA payments in IMPLAN plus NBER TAXSIM rates
parameter tl0(r) tax rate - labor /
$ondelim
$include 'build\data\income_tax_rates.csv'
$offdelim
/;
tl0(r) = tl0(r)+sum(pub, fs(r,pub,"sstw","empl"))/sum(h, le0(r,h))
               +sum(pub, fs(r,pub,"sstf","empl"))/sum(h, le0(r,h));

* capital tax rate
* from EPPA (Gurgel et al., 2007)
tk0(r) = 0.376;

* consumption tax rate - read from csv file
parameter tc0(r) tax rate - consumption /
$ondelim
$include 'build\data\consumption_tax_rates.csv'
$offdelim
/;

* rescale labor demand and endowments to market prices based on tax rate
ld0(r,s) = ld0(r,s)/(1+tl0(r));
le0(r,h) = le0(r,h)/(1+tl0(r));

* rescale capital demand and endowments to market prices based on tax rate
kd0(r,s) = kd0(r,s)/(1+tk0(r));
re0(r,h) = re0(r,h)/(1+tk0(r));



* ------------------------------------------------------------------------------
* disaggregate the oil and gas extraction data
* ------------------------------------------------------------------------------

* the initial aggregation in "build/aggregation_map/default_aggregation.gms"
* includes oil and gas extraction in the crude oil sector cru and only natural
* gas distribution in the sector gas. this a result of the minimum aggregation
* level in the underlying IMPLAN data. therefore, this script disaggregates the
* oil and gas extraction in cru into natural gas and crude oil extraction,
* shifting the natural gas extraction to the gas sector. to determine natural
* gas's share of consumption/use we assume that crude oil serves as an
* intermediate input only to the petroleum refining sector and intermediate
* inputs to all other sectors represent natural gas, as does household and
* government consumption and investment demand. some of the intermediate input
* to the refining sector will also be natural gas. to determine that share and
* natural gas's share of production and trade we minimize the sum of squared
* deviations for those shares from observed values or assumed shares conditional
* on market clearance conditions and the assumption of weakly positive domestic
* use of production. The observed or assumed shares we try to match are derived
* as follows:
*
* 1.) the observed share of natural gas production by region is defined using
*     EIA data on crude oil and natural gas production by state aggregated up to
*     the regional level. to arrive at a value share we multiply state level
*     production quantities by EIA data on state level wellhead prices for crude
*     oil and city gate natural gas prices as a proxy for natural gas well head
*     prices.
*
* 2.) the share of natural gas international imports and exports by region are
*     defined using census data on state level international imports and exports
*     of crude oil and natural gas aggregated to the regional level.
*
* 3.) a region's intrantional import share of natural gas is assumed to similar
*     to the region's share of natural gas use relative to the crude oil and
*     natural gas total. a region's intranational export share of natural
*     gas is assumed to be similar to the share of natural gas production in the
*     region.
*
* 4.) the observed share of natural gas used as an intermediate input in the
*     refining sector is estimated based on national annual averages of crude
*     oil and natural gas inputs to the sector collected by EIA and converted
*     to values using the Brent and Henery Hub average annual prices as reported
*     by EIA.
*
* exogenous data on crude oil and natural gas extraction are generated by
* "build/get_oil_and_gas_data.R".
*

sets
  s_gas(s)              natural gas sectors                             / 21 /
  s_cru(s)              natural gas and crude oil extraction sectors    / 20 /
  s_ref(s)              petroleum refining sector                       / 156 /;

parameters
  y_og(r,*)             observed - natural gas and crude oil production
  x_og(r,*)             observed - natural gas and crude oil international exports
  m_og(r,*)             observed - natural gas and crude oil international imports
  obs_gas_y_share(r)    observed - natural gas share of production
  obs_gas_x_share(r,trd) observed - natural gas share of international exports
  obs_gas_m_share(r,trd) observed - natural gas share of international imports
  obs_ref_gas_share     observed - natural gas share of refining sector input;

variables
  gas_y_share(r)        benchmark - natural gas share of production
  gas_x_share(r,trd)    benchmark - natural gas share of exports
  gas_m_share(r,trd)    benchmark - natural gas share of imports
  ref_gas_share(r)      benchmark - natural gas share of refining sector input
  ssd                   objective - sum of squared deviations
  use_gas(r)            regional use of natural gas;

equations
  eq_use_gas(r)         regional use of natural gas
  mc_gas(r)             market clearance - natural gas
  mc_dtrd               market clearance - domestic trade
  dm_gas(r)             (weakly) positive domestic own use
  dm_cru(r)             (weakly) positive domestic own use
  obj                   objective - sum of squared deviations;


* load the additional data on oil and gas production and trade
$include "build/data/oil_and_gas_data.gms"

* observed - natural gas share of production
obs_gas_y_share(r)$(y_og(r,"gas")+y_og(r,"cru"))
  = y_og(r,"gas")/(y_og(r,"gas")+y_og(r,"cru"));

* observed natural gas share of refining sector input
* based on eia data
obs_ref_gas_share = (3*881611*1000)/(43.29*5924395*1000);

* observed - natural gas share of international imports
obs_gas_m_share(r,"ftrd")$(m_og(r,"gas")+m_og(r,"cru"))
  = m_og(r,"gas")/(m_og(r,"gas")+m_og(r,"cru"));

* observed - natural gas share of international exports
obs_gas_x_share(r,"ftrd")$(x_og(r,"gas")+x_og(r,"cru"))
  = x_og(r,"gas")/(x_og(r,"gas")+x_og(r,"cru"));

* observed - natural gas share of intranational imports
obs_gas_m_share(r,"dtrd") = 1-(1-obs_ref_gas_share)*sum((s_cru,s_ref), id0(r,s_cru,s_ref))
                               /sum(s_cru, sum(s, id0(r,s_cru,s))+sum(h, cd0(r,s_cru,h))
                                 +g0(r,s_cru)+i0(r,s_cru));

* observed - natural gas share of intranational exports
obs_gas_x_share(r,"dtrd") = obs_gas_y_share(r);

* regional use of natural gas
eq_use_gas(r)..
  use_gas(r) =e= sum(s_cru,  sum(s, id0(r,s_cru,s))
                            -(1-ref_gas_share(r))*sum(s_ref, id0(r,s_cru,s_ref))
                            +sum(h, cd0(r,s_cru,h))
                            +g0(r,s_cru)
                            +i0(r,s_cru));

* market clearance - natural gas
mc_gas(r)..
  sum(s_cru, sum(trd, gas_m_share(r,trd)*m0(r,s_cru,trd))
             +gas_y_share(r)*y0(r,s_cru))
    =e= use_gas(r)+sum(s_cru, sum(trd, gas_x_share(r,trd)*x0(r,s_cru,trd)));

* market clearance - domestic trade
mc_dtrd..
  sum(r, gas_m_share(r,"dtrd")*sum(s_cru, m0(r,s_cru,"dtrd"))
        -gas_x_share(r,"dtrd")*sum(s_cru, x0(r,s_cru,"dtrd"))) =e= 0;

* (weakly) positive domestic own use
dm_gas(r)..
  .98*gas_y_share(r)*sum(s_cru, y0(r,s_cru))
    =g= sum(trd, gas_x_share(r,trd)*sum(s_cru, x0(r,s_cru,trd)));

* (weakly) positive domestic own use
dm_cru(r)..
  .98*(1-gas_y_share(r))*sum(s_cru, y0(r,s_cru))
    =g= sum(trd, (1-gas_x_share(r,trd))*sum(s_cru, x0(r,s_cru,trd)));

* objective - sum of squared deviations
obj..
  ssd =e= sum((r,trd), power(gas_x_share(r,trd)-obs_gas_x_share(r,trd),2)
                      +power(gas_m_share(r,trd)-obs_gas_m_share(r,trd),2)
                      +power(gas_y_share(r)-obs_gas_y_share(r),2)
                      +power(ref_gas_share(r)-obs_ref_gas_share,2));

* lower bound for shares
gas_m_share.lo(r,trd) = 0;
gas_x_share.lo(r,trd) = 0;
ref_gas_share.lo(r)   = 0;

* upper bound for shares
gas_m_share.up(r,trd) = 1;
gas_x_share.up(r,trd) = 1;
ref_gas_share.up(r)   = 1;

* adjust natural gas production and trade shares -> there are no oil or gas
* wells in some regions and therefore no production data available for
* either commodity so the included production share is zero. the non-zero output
* value in the sam for this sector is likely due to services that are exported
* or processing of imported commodities. since natural gas accounts for nearly
* all foreign exports and regional use in these regions in the default
* aggregation (i.e., new england) we assume that all of the regional production
* value is associated with natural gas. crude oil exports in the data are
* deminimus but to match this production assumption we zero them out completely.
gas_y_share.fx(r)$(obs_gas_y_share(r) eq 0) = 1;
gas_x_share.fx(r,"ftrd")$(obs_gas_y_share(r) eq 0) = 1;

* define the disaggregation problem
model cal_gas / eq_use_gas mc_gas mc_dtrd dm_gas dm_cru obj /;

* solve the disaggregation problem
solve cal_gas using nlp minimizing ssd;

loop(s_cru,

* disaggregate trade
m0(r,s_gas,trd) = m0(r,s_gas,trd)+gas_m_share.l(r,trd)*m0(r,s_cru,trd);
x0(r,s_gas,trd) = x0(r,s_gas,trd)+gas_x_share.l(r,trd)*x0(r,s_cru,trd);
m0(r,s_cru,trd) = (1-gas_m_share.l(r,trd))*m0(r,s_cru,trd);
x0(r,s_cru,trd) = (1-gas_x_share.l(r,trd))*x0(r,s_cru,trd);

* partition the input into refining into a portion that comes through gas
* and have the remaining portion come through cru
id0(r,s_gas,s_ref) = id0(r,s_gas,s_ref)+ref_gas_share.l(r)*id0(r,s_cru,s_ref);
id0(r,s_cru,s_ref) = (1-ref_gas_share.l(r))*id0(r,s_cru,s_ref);

* for non-refining sectors crude oil inputs are assumed to be zero so all
* inputs should come through gas
id0(r,s_gas,s)$(not s_ref(s)) = id0(r,s_gas,s)+id0(r,s_cru,s);
id0(r,s_cru,s)$(not s_ref(s)) = 0;

* all investment demand is in the form of natural gas - this matters
i0(r,s_gas) = i0(r,s_gas)+i0(r,s_cru);
i0(r,s_cru) = 0;

* all government demand is in the form of natural gas - this shouldn't be
* necessary but kept in for completeness
g0(r,s_gas) = g0(r,s_gas)+g0(r,s_cru);
g0(r,s_cru) = 0;

* all consumption demand is in the form of natural gas - this shouldn't be
* necessary but kept in for completeness
cd0(r,s_gas,h) = cd0(r,s_gas,h)+cd0(r,s_cru,h);
cd0(r,s_cru,h) = 0;

* partition intermediate inputs based on share of natural gas production
id0(r,s,s_gas) = id0(r,s,s_gas)+gas_y_share.l(r)*id0(r,s,s_cru);
id0(r,s,s_cru) = (1-gas_y_share.l(r))*id0(r,s,s_cru);

* partition capital demand based on share of natural gas production
kd0(r,s_gas) = kd0(r,s_gas)+gas_y_share.l(r)*kd0(r,s_cru);
kd0(r,s_cru) = (1-gas_y_share.l(r))*kd0(r,s_cru);

* partition labor demand based on share of natural gas production
ld0(r,s_gas) = ld0(r,s_gas)+gas_y_share.l(r)*ld0(r,s_cru);
ld0(r,s_cru) = (1-gas_y_share.l(r))*ld0(r,s_cru);

);

* production tax rate for the natural gas sector needs to be adjusted to reflect
* the extraction activities that would be taxed at that rate - done by keeping
* solving for the tax rate that keeps output value equivalent to the original
ty0(r,s_gas) = (1-(sum(ss, id0(r,ss,s_gas))+(1+tl0(r))*ld0(r,s_gas)
                  +(1+tk0(r))*kd0(r,s_gas))
                 /(y0(r,s_gas)+gas_y_share.l(r)*sum(s_cru, y0(r,s_cru)))
                  )$(y0(r,s_gas)+gas_y_share.l(r)*sum(s_cru, y0(r,s_cru)));



* ------------------------------------------------------------------------------
* load the benchmark household data
* ------------------------------------------------------------------------------

$include "build/data/population_data.gms"



* ------------------------------------------------------------------------------
* check balance of national sam
* ------------------------------------------------------------------------------

set
  sam_accounts          macro sam accounts / activity, commodity, labor,
                                             capital, household, government,
                                             investment, national, foreign /;

alias (sam_accounts,aa);

* benchmark - government to household lump sum transfer
* set to the level that balances the household budget constraint
incadj0(r,h) = sum(s, (1+tc0(r))*cd0(r,s,h))+invh0(r,h)-re0(r,h)-le0(r,h)
               -bopdef0(r,h);

parameters
  macro_sam             macro sam;

* build macro sam
macro_sam("activity",   "commodity"  ) = sum((r,s), y0(r,s));
macro_sam("commodity",  "activity"   ) = sum((r,s,ss), id0(r,ss,s));
macro_sam("commodity",  "household"  ) = sum((r,s,h), cd0(r,s,h));
macro_sam("commodity",  "government" ) = sum((r,s), g0(r,s));
macro_sam("commodity",  "investment" ) = sum((r,s), i0(r,s));
macro_sam("commodity",  "national"   ) = sum((r,s), x0(r,s,"dtrd"));
macro_sam("commodity",  "foreign"    ) = sum((r,s), x0(r,s,"ftrd"));
macro_sam("labor",      "activity"   ) = sum((r,s), ld0(r,s));
macro_sam("capital",    "activity"   ) = sum((r,s), kd0(r,s));
macro_sam("household",  "labor"      ) = sum((r,h), le0(r,h));
macro_sam("household",  "capital"    ) = sum((r,h), re0(r,h));
macro_sam("household",  "government" ) = sum((r,h), incadj0(r,h));
macro_sam("household",  "foreign"    ) = sum((r,h), bopdef0(r,h));
macro_sam("government", "activity"   ) = sum((r,s),  ty0(r,s)*y0(r,s)
                                                    +tk0(r)*kd0(r,s)
                                                    +tl0(r)*ld0(r,s));
macro_sam("government", "household"  ) = sum((r,s,h), tc0(r)*cd0(r,s,h));
macro_sam("investment", "household"  ) = sum((r,h), invh0(r,h));
macro_sam("national",   "commodity"  ) = sum((r,s), m0(r,s,"dtrd"));
macro_sam("foreign",    "commodity"  ) = sum((r,s), m0(r,s,"ftrd"));

* compute the row and colum totals and check the differences
macro_sam(sam_accounts,"total") = sum(aa, macro_sam(sam_accounts,aa));
macro_sam("total",sam_accounts) = sum(aa, macro_sam(aa,sam_accounts));
macro_sam("r-c diff",sam_accounts) =  macro_sam(sam_accounts,"total")
                                     -macro_sam("total",sam_accounts);

display macro_sam;

parameters
  macro_stats(*,*)      basic statistics for the macro economy;

macro_stats("consumption exp","value") = sum((r,h,s), cd0(r,s,h));
macro_stats("investment exp","value")  = sum((r,s),i0(r,s));
macro_stats("government exp","value")  = sum((r,s), g0(r,s));
macro_stats("exports","value")         = sum((r,s), x0(r,s,"ftrd"));
macro_stats("imports","value")         = sum((r,s), m0(r,s,"ftrd"));
macro_stats("gdp by exp","value")      =  macro_stats("consumption exp","value")
                                         +macro_stats("investment exp","value")
                                         +macro_stats("government exp","value")
                                         +macro_stats("exports","value")
                                         -macro_stats("imports","value");
macro_stats("labor va","value")        = sum((r,s), (1+tl0(r))*ld0(r,s));
macro_stats("capital va","value")      = sum((r,s), (1+tk0(r))*kd0(r,s));
macro_stats("prod tax","value")        = sum((r,s), ty0(r,s)*y0(r,s));
macro_stats("gdp by va","value")       =  macro_stats("labor va","value")
                                         +macro_stats("capital va","value")
                                         +macro_stats("prod tax","value");
macro_stats("households","value")      = sum((r,h), n0(r,h));

display macro_stats;



* ------------------------------------------------------------------------------
* save benchmark dataset
* ------------------------------------------------------------------------------

execute_unload 'build\data\noaggr_benchmark.gdx', f,s,r,h,t,id0,ld0,kd0,cd0,x0,
                                                  m0,i0,g0,re0,le0,invh0,n0
                                                  ty0,tl0,tc0,tk0,trd;
